% Example,
%   % Read two greyscale images of Lena
clear all;
close all;
% Imoving=imread('images/modtest1.png');
% Istatic=imread('images/modtest1-15d.png');
Imoving=imread('images/lenag1large.png');
Istatic=imread('images/lenag2large-15d.png');

%
%   % Register the images
options = struct('Similarity','p','Registration','NonRigid','SigmaFluid',4,...
    'Alpha',4,'SigmaDiff',1,'Interpolation','Linear','maxiter',500);
[Ireg,Bx,By,Fx,Fy] = affineregister_images(Imoving,Istatic,1,1,options);
%[Ireg,Bx,By,Fx,Fy] = affineregister_images(Imoving,Istatic,struct('Similarity','p','Registration','NonRigid'));
%
%   % Show the registration result
%imwrite(Ireg,'result_DemonswithAF.png');
%imwrite(Ireg,'result_AFbuildin2.png');
%imwrite(Ireg,'result_standardDemonsbox.png');
%[X,Y] = ndgrid(linspace(0,1,size(Imoving,1)),linspace(0,1,size(Imoving,2)));
[X,Y]=ndgrid(0:(size(Imoving,1)-1),0:(size(Imoving,2)-1));

figure,
% subplot(2,2,1), imshow(Imoving); title('moving image');
% subplot(2,2,2), imshow(Istatic); title('static image');
% subplot(2,2,3), imshow(Ireg); title('registerd moving image');
subplot(1,3,1), imshow(Imoving); title('moving image');
subplot(1,3,2), imshow(Istatic); title('static image');
subplot(1,3,3), imshow(Ireg); title('registerd moving image');

hold on
contour( X+Bx, 10, 'color', 'r' );
contour( Y+By, 10, 'color', 'r' );
hold off

% [A, b] = GetAffineMatrix(X+Bx, Y+By, X, Y);
% Q = [A(1,1) A(1,2) b(1);
%      A(2,1) A(2,2) b(2);
%      0      0      1];
% %   % Show also the static image transformed to the moving image
% Ireg2=movepixels(Istatic,Fx,Fy,[], 3);
% 
% subplot(2,2,4), imshow(Ireg2); title('registerd static image');
%
%  % Show the transformation fields
figure,
subplot(2,2,1), imshow(Bx,[]); title('Backward Transf. in x direction');
subplot(2,2,2), imshow(Fx,[]); title('Forward Transf. in x direction');
subplot(2,2,3), imshow(By,[]); title('Backward Transf. in y direction');
subplot(2,2,4), imshow(Fy,[]); title('Forward Transf. in y direction');
%
% % Calculate strain tensors
E = strain(Fx,Fy);
% % Show the strain tensors
% figure,
% subplot(2,2,1), imshow(E(:,:,1,1),[]); title('Strain Tensors Exx');
% subplot(2,2,2), imshow(E(:,:,1,2),[]); title('Strain Tensors Exy');
% subplot(2,2,3), imshow(E(:,:,2,1),[]); title('Strain Tensors Eyx');
% subplot(2,2,4), imshow(E(:,:,2,2),[]); title('Strain Tensors Eyy');

disp(['SAD = ' num2str(sum(sum(abs(Istatic-Ireg))))]);
mse = sum(sum((Istatic-Ireg).^2))/(size(Istatic,1)*size(Istatic,2));
disp(['MSE = ' num2str(mse)]);

dIm1 = Istatic - mean(Istatic(:));
dIm2 = Ireg - mean(Ireg(:));
numerator = sum(sum(dIm1.*dIm2));
denominator = sqrt(sum(sum(dIm1.^2)))*sqrt(sum(sum(dIm2.^2)));

ncc = numerator/denominator;
disp(['NCC = ' num2str(ncc)]);